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Linear augmentation has recently been shown to be effective in targeting desired stationary solu¬ 
tions, suppressing bistablity, in regulating the dynamics of drive response systems and in controlling 
the dynamics of hidden attractors. The simplicity of the procedure is the highlight of this scheme 
but at the same time questions related to its general applicability still need to be addressed. Fo¬ 
cusing on the issue of targeting stationary solutions, this work demonstrates instances where the 
scheme fails to stabilize the required solutions and leads to other complicated dynamical scenarios. 
Appropriate examples from conservative as well as dissipative systems are presented in this regard 
and potential applications for relevant observations in dissipative predator-prey systems are also 
discussed. 
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I. INTRODUCTION 

Studies on coupled nonlinear systems have explored a 
wide variety of emergent dynamical phenomena, namely 
synchronization [I|, oscillator suppression J3L mul- 
tistability 0, hysteresis extreme-events [313 
which can be exploited in applications, to model nat¬ 
ural phenomena or in regulating the system behavior 
for instance. Controlling dynamical systems towards a 
desired behavior is an important research topic in non¬ 
linear sciences Q- Starting with chaos control &li3> 
research in this domain now also extends towards con¬ 
trol of multistability [l^, patterns and spatio-temporal 
chaos BUI, noisy systems [I1,[I3 , methods of sta¬ 
bilizing unstable stationary solutions [3,1 [13 etc. A 
greater understanding of these different regulatory as¬ 
pects have greatly contributed towards development of 
related novel and highly efficient procedures. Consider¬ 
ing noninvasive (without changing the intrinsic system 
parameters) mechanisms leading to stabilization of sta¬ 
tionary solutions, oscillator suppression via coupling non¬ 
linear systems has been discussed extensively in literature 
(see Refs. for detailed reviews). This suppression 

is majorly observed as a consequence of parameter het¬ 
erogeneity between the coupled units [l9l - [^. p resence 
of time-delayed [13, [13 /conjugate variables 123 in the 
coupling function or through dynamic coupling |25l |. 

Recently, linear augmentation has also been suggested 
as another practical alternative leading to oscillator sup¬ 
pression, achieved by coupling systems to a linear feed¬ 
back consisting of a simple decaying function [13 . Lately, 
studies have also used linear augmentation effectively for 
controlling bistability [13, dynamics of a drive response 
system |28l | and in controlling the dynamics of hidden 
attractors [13- With respect to stabilizing stationary 
solutions. Ref. [13 discussed results for an augmented 
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Lorenz system where either the stationary solutions of 
the original system or those of the augmented system 
could be stabilized by picking an appropriate feedback 
function; former being quite relevant from an applica¬ 
tion perspective. The paper also presented some param¬ 
eter space scans highlighting the regimes where linear 
augmentation works and where it does not, which al¬ 
though is instructive but is also very system specific at 
the same time. At this point, one must question the abil¬ 
ity of linear augmentation towards stabilizing the sta¬ 
tionary solutions in a more general sense, namely the 
systems, parameter settings and coupling configurations 
where the scheme works and where it does not? In this 
paper, we will look at some simple examples of linearly 
augmented systems demonstrating the fact that we need 
to be quite careful before picking linear augmentation in 
applications. These examples illustrate that there could 
be situations where even picking an appropriate feedback 
function does not guarantee that the required stationary 
solutions will be necessarily stabilized. We will see that 
the mechanism appears to be highly dependent on the in¬ 
trinsic properties of the oscillators in consideration, the 
stationary solutions we want to target, and also on how 
these systems are augmented by/coupled to the feedback. 
Furthermore, we will also discuss instances where the fail¬ 
ures of the procedure can be exploited in applications. 

The manuscript is arranged as follows: Linear augmen¬ 
tation is introduced in the following Sec.llTl In Sec. lIIIl we 
will have a look at results for linearly augmented conser¬ 
vative and dissipative dynamical systems. In Sec. IIII Al 
results for partially and fully augmented harmonic oscil¬ 
lator are presented. Sec. ME discusses results for con¬ 
servative Duffing oscillator under similar augmentations, 
and in Sec IIII Cl we will have a look at the behavior of 
augmented dissipative population models and also briefly 
discuss possible applications for certain observations in 
these systems. The manuscript concludes with a sum¬ 
mary of results in Sec. IlYl Additional details on certain 
dynamical aspects of harmonic oscillator and Duffing sys¬ 
tem which were excluded from main manuscript are pro- 









2 


vided in Appendices 0 [B1 and [Cl 


II. LINEAR AUGMENTATION 

General representation of a linearly augmented dynam¬ 
ical system is, 


X = f(x) -I- eu 
u = —ku — e(x — b)^. 

where the column vector x = [a;i, 0:2,..., S 

([...]^ corresponds to the transpose) contains the sys¬ 
tems variables, and u is the augmentation variable, e = 
[ei, £2, ■ ■ ■, G is a column vector with infor¬ 

mation regarding the coupling strength of the interac¬ 
tion between the dynamical variables and u; augmenta¬ 
tion/coupling term corresponding to the component 
Xi is £iU V i = 1 , 2 ,..., N, and = 0 if Xi is not coupled 
to u. b £ is an arbitrary vector and k is the decay 
constant [s^ which could be used to control the transient 
time leading to stabilization of stationary solutions (^ . 
Vector b = X* = [xi*,X2*, ■ ■ ■ ,xn*]^ £ where x* 
satisfies x|x=x* = f(x*) = 0 if we want to stabilize a sta¬ 
tionary solution X* of the original system. Substituting 
a value of b ^ x* can stabilize stationary solutions of 
augmented system for which X = [xi,X2, ■ ■ ■, x'n, £ 
= 0. The term £(x — b)^ gives the dot product of 
the corresponding column vectors. 

In the following, we will look at examples of augmented 
conservative as well as dissipative dynamical systems 
which highlight the limitations of this procedure. The 
terms augmentation/augmented and coupling/coupled 
have been used synonymously in the following text. 



III. EXAMPLES 

Here we will discuss the instances of systems controlled 
via linear augmentation. We will first look at two exam¬ 
ples of conservative systems, namely the harmonic os¬ 
cillator and conservative Duffing oscillator where linear 
augmentation will be used to stabilize their stationary 
solutions. 


A. Harmonic oscillator 

Equations describing a linearly augmented harmonic 
oscillator are: 


X = y + eu, 
y = -uj'^x, 

ii = —ku — ex, (2) 

where lj is the frequency of the oscillator, k is the aug¬ 
mentation parameter, and e is the coupling strength. The 
first two equations governing the evolution of x and y 



FIG. 1. (Color online) a) Bifurcation diagram (black dots), 
largest eigenvalues (red symbols) and the time series of x 
variable before and after the stabilization of origin as inset. 
Corresponding phase space plots are in b) (no augmentation: 
£ = 0) and c) (with augmentation: e = 0.5). d) Variation in 
the eigenvalues for higher values of augmentation strength e. 
Parameter values for calculation are lj = 2 and fc = 2. 


correspond to the original harmonic oscillator dynamics. 
In the absence of augmentation, harmonic oscillator con¬ 
serves total energy, which stays constant on the ellipses 
shown in Fig. |l]b). Each of these ellipses correspond to 
the system evolution following different initial values of 
position and momentum, and hence, different conserved 
total energies. Harmonic oscillator has a;* = 0, y* = 0 
as the only stationary solution and note that the aug¬ 
mentation term only appears in the rate equation of the 
position variable x at this point. In case of a successful 
stabilization, the required stationary solution of the full 
system should be {x*,y*,u*) = ( 0 , 0 , 0 ) (origin) where 
the system effectively decouples from augmentation. 

The characteristic eigenvalue equation at the origin is, 

(A + k)(X^+uj^)+e^A = 0. (3) 

For £ = 0, the eigenvalues for the full system are Ai ^2 = 
zLicu which correspond to the non-hyperbolic stationary 
solution at the origin, and A 3 = — fc corresponding to the 
decay of the augmentation variable u; which evolves as 
u(t) oc exp(—kt) in this case. The parameter values for 
the following calculations were fixed at w = 2, and fc = 2. 
For the evolution of the augmented system (e > 0), the 
bifurcation diagram of the system with increasing e 
values is shown in Fig. [T]a) (black dots). It is seen that 
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with an increasing e, the system which was conservative 
for e = 0 becomes dissipative and gets into a stable origin 
regime even for quite small values of e. Rewriting Eq. [3] 
as, 


+ kX^ + + e^)A + kuj^ = 0, 


(4) 


and applying the Routh-Hurwitz criteria (RHC) [32|, we 
can deduce that the roots of this equation are all either 
negative or have negative real parts (in case of complex 
roots) V e > 0. Largest eigenvalues of the Jacobian (red 
symbols) are also plotted along with the bifurcation dia¬ 
gram in Fig. [T]a) which demonstrate the transition from 
oscillatory to stationary state for e > 0. Considering 
the behavior of this system for large e, we can see that 
the largest eigenvalue A = Ai —>■ 0 from Eq. S] in this 


limit. Since the discriminant 33| of the cubic character¬ 
istic Eq.|4]is negative V e > 0, this implies that the system 
has one real eigenvalue and a pair of complex conjugate 
eigenvalues in this range. The negative real part of these 
complex eigenvalues for large e can therefore be estimated 
by equating the sum of all eigenvalues to the trace of the 
Jacobian tr(J), giving Re(A 2 , 3 ) = —A:/2. This further 
implies that as e —^ oo, convergence to the origin gets 
slower although origin is stable in the entire £ > 0 range 
and any change in stability will only occur as £ —>• oo 
when Ai = 0. 

Now let us consider a more general case of an aug¬ 
mented harmonic oscillator given by, 


X = y + £iM, 

y =+ e 2 U, (5) 

ii = —ku — Six — £ 22 /, 

where the augmentation now appears in the rate equa¬ 
tions of both position x and momentum y with coupling 
strengths £i, £2 respectively. The eigenvalue equation in 
this case is, 

(A -f k){X^ w^) -|- A(£i^ -|- £ 2 ^) + £i£ 2(1 ~ w^) = 0. 

( 6 ) 

Substituting £2 = 0 and £1 = £ in Eq. [S] yields the dy¬ 
namics of Eq. [21 Similarly, for £1 = 0 and £2 = £, we 
obtain a case where the system is only coupled in the y 
variable for which the characteristic Eq.jSlis exactly iden¬ 
tical to Eq. 01 and therefore the stability characteristics 
of the origin are identical and independent of whether the 
system is augmented m x or y. In the previous example, 
we saw a situation where linear augmentation success¬ 
fully stabilized the origin for the entire range of £ > 0. 
Now considering £1 = £2 = £ (the system is similarly 
augmented in both variables), in which case Eq. [3] gives, 

A^ -I- kX^ -I- (w^ -I- 2£^)A -I- -I- £^(1 — w^) = 0. (7) 

Using the RHC, it can be seen that this equation will have 
all negative eigenvalues iff -I- £^(1 — w^) > 0 which 
gives a stability regime of0<£<£*Va;>l where 

I 

£* = u}\ — -, and for higher values of £, RHC suggests 

Vo;"' — ! 



FIG. 2. (Color online) Largest eigenvalue estimates for sta¬ 
tionary solutions {x*,y*) = (0,0) (black), {x*,y*) = (±1,0) 
(red) for partially augmented Duffing system. 


appearance of positive eigenvalue/eigenvalues. For large 
£, we can get an estimate of largest eigenvalue Ai —>■ 

^ > 0 V w > 1. Since the discriminant is negative 


V £ > 0, the remaining complex conjugate eigenvalue pair 
have a negative real part given by Re(A 2 , 3 ) = (Tr(J) — 

Ai)/2 = — ± ^ ^ j. Therefore, unlike in the 


previous example, we can see that origin here is unstable 
for large e. The expression for e* also shows that a higher 
value of k can extend the coupling range for a stable 
origin. This result is the first instance of unexpected 
behavior as we would normally expect a higher coupling 
value to keep the origin stable. Furthermore, in the e > 
e* regime it is numerically observed that the trajectories 
escape to infinity which is also quite unexpected. 


One of the primary reasons behind considering aug¬ 
mented harmonic oscillator in this study is the fact that 
it is highly solvable and therefore can provide necessary 
insights into the physical mechanisms behind the desir¬ 
able as well as undesirable behaviors. It turns out that in 
case of a successful stabilization, this system represents 
a forced harmonic oscillator where the steady state so¬ 
lution (which is completely determined by the forcing) 
decays to the origin along with an exponentially decay¬ 
ing force. Similarly, the case where the trajectories es¬ 
cape to infinity corresponds again to a forced system but 
this time being driven by an exponentially diverging force 
which is analogous to a situation where energy is being 
pumped into the system. Therefore the steady state so¬ 
lution in this case diverges along with the diverging force 
explaining the unexpected behavior of escaping trajecto¬ 
ries observed for the fully augmented system. Details of 
the calculations leading to these deductions are available 
in Appendix El 
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B. Duffing oscillator 

General equations for a linearly augmented Duffing os¬ 
cillator with no damping or forcing can be written as: 

X = y + eiu, 
y = X — + e 2 U, 

u = —ku — ei{x — X*) — £ 2 ( 2 / — y*). ( 8 ) 

Uncoupled Duffing system has an invariant of motion 
(also the Hamiltonian of the system) given by H{x, y) = 
y^ 12 — x^ 12 -b x'^/4 and stationary solutions: {x*,y*) = 
(±1,0), (0,0). The trajectories of this system evolve on 
the double well potential surface of H{x, y) starting from 
different initial conditions for £1 = £2 = 0. Siirrilar to 
the previous example, for a successful stabilization, the 
required stationary solutions of the full system should 
be {x*,y*,u*) = ( 0 , 0 , 0 ) or (± 1 , 0 , 0 ) where the system 
effectively decouples from the augmentation. These so¬ 
lutions will be referred to as {x*,y*) = ( 0 , 0 ) (origin) or 
(± 1 , 0 ) in the following. 

For the system in Eq. [51 the characteristic eigenvalue 
equation can be expressed as, 

(A ± k)[X^ ± 3a;* — 1) ± A(£i^ ± £ 2 ^) 

±£i£ 2 ( 2 - 3 a;*^) =0. (9) 

Now similar to the harmonic oscillator example, consid¬ 
ering partial augmentation with £ 1 ( 2 ) = £, and £ 2 ( 1 ) = 0 
first, Eq. [5] suggests that the stability characteristics for 
the stationary solutions are again independent of whether 
the system is being augmented in x or y. For this partial 
augmentation, Eq. [9] gives, 

{X + k){X^+3x*^-1)+Xe^ = 0. (10) 

Substituting x* = 0,±1, and rearranging the terms, we 
can obtain the characteristic eigenvalue equations for 
these stationary solutions as, 

X^+ kX^+ {e^-l)X-k = 0, (11) 

for {x*,y*) = ( 0 , 0 ) (hyperbolic for £ = 0 ), and 

X^+ kX‘^+ {e^+ 2)X + 2k = 0, (12) 

for {x*,y*) = (± 1 , 0 ) (non hyperbolic for £ = 0 ) respec¬ 
tively. It is straightforward to check that the largest 
eigenvalue Ai — 0 for larger e values in both these cases 
which implies that a stable/unstable stationary solution 
will retain its stability characteristics until a stability 
change (zero crossing of the eigenvalue/s) occurs in the 
£ —> 00 limit. Furthermore using the RHC, it is easily 
verifiable that Eq.[TT]will always have positive root/roots, 
whereas Eq. [T^] will have all negative roots V £ > 0; 
which implies that the a:* = 0 is always unstable and 
a:* = ±1 is always stable. Therefore, we see that partial 
augmentation works for stabilizing {x*,y*) = (± 1 , 0 ) but 
fails completely to stabilize the origin {x*,y*) = ( 0 , 0 ). 
Fig. m shows the largest eigenvalue calculations which 
verify these deductions. This brings us to an important 



FIG. 3. (Color online) Figs, a), b) and c) show the bifur¬ 
cation diagrams (black dots) and the largest eigenvalues (red 
symbols) for {x*,y*) = (1,0), (0,0), and (—1,0) respectively. 
Related Figs, a.l: for e = 0.4, the system is bistable and 
the two related transient behaviors (in blue and green and 
likewise for other cases), a.2: for e = 1, the trajectory ap¬ 
proaching the stable stationary solution (1, 0), and a.3 shows 
an arbitrary time series for e = 2.5. Similarly in b.l: bistabil¬ 
ity, and in b.2: the system approaching the stable stationary 
solution (0,0) is shown. Identically, c.l, c.2, and c.3 show the 
bistability (e = 0.4), stabilization of (—1,0) (e = 1) and an 
arbitrary time series at e = 2.5. 

observation that there might exist situations where it is 
not possible to target the required stationary solution 
even on using an appropriate feedback function with any 
combination of k and e values. 

Now considering identical augmentation with £1 = 
£2 = £ and we will see that this system has some in- 
teresting properties. Fig. |5| shows the bifurcation dia¬ 
grams [3l| of the system as we try targeting different 
stationary solutions: For {x*,y*) = (1,0), the bifurca¬ 
tion diagram (black dots) is shown in Fig. |5|a). Appro¬ 
priate transient trajectories in different coupling regimes 
are also shown in related Figs. |3| (a.l), (a.2), (a.3). It 
is observed that even for very small coupling values, the 
system quickly gets into a stable stationary state regime, 
although, for smaller values of £, it exhibits bistability. 
The transient trajectories in this parameter regime are 
shown in Fig. |5| (a.l). We observe that the augmen¬ 
tation is stabilizing our desired stationary solution at 
{x*,y*) = (1,0), but along with it, other stationary so¬ 
lutions which are £ dependent are also getting stabilized 
on starting with different initial conditions. These other 
stationary solutions for the augmented system here are 
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given by, 



y ± = X ±-X ± 

* y*± 

u ± = -, 


(13) 


and solutions {x*-,y*z*-) are observed to coexist 
along with {x*,y*) = (1,0). For higher coupling val¬ 
ues, bistability terminates via a saddle node bifurca¬ 
tion when the stable branch of stationary solutions 
{x*-,y*_,u*-) collides with the unstable branch of 
(a;* + , 2 /* + , M*+) (circles) as shown in in Fig. [3] (a) at 
[k 

SSN = \ The system also exhibits hysteresis in this 

bistable regime and a brief discussion regarding this ob¬ 
servation is available in Appendix [BJ Beyond this regime 
for a range of values in e > esN, (x* ,y*) = (1,0) remains 
as the only stable attractor as shown in Fig. [3] (a),(a.2). 

In absence of augmentation, the eigenvalues for 
{x*,y*) = (1,0) are complex: Ai ^2 = ±*-\/2. For the aug¬ 
mented system, the characteristic equation can therefore 
be written as: 


X^ + kX^ + 2X{l+e^) + {2k-e^)=0. (14) 

The RHC shows that this equation will have all negative 
roots for {2k —e^) > 0 and positive root/roots appear for 
e > This gives us the transition threshold for the 

destabilization of the stationary solution as e* = 
at which the eigenvalue/s cross the zero axis. Since the 
discriminant is negative, the characteristic equation has 
one real and two complex conjugate roots. Considering 
large e behavior, it is seen that the largest eigenvalue 
Ai —>■ 1/2 which implies that {x*,y*) = (1,0) is unstable 
in this range. The real part of the remaining complex 
conjugate eigenvalue pair is Re(A 2 , 3 ) = —{2k + l)/4. At 
e*, Eq. [T4]can be rewritten as, 

X{X‘^+ kX + 2{l + 2k)) = 0, (15) 

which consequently gives the eigenvalues as Ai = 0 and 
A2,3 = {—k ± — 8(1 -I- 2k))/2. We can see that A2,3 

will be a complex conjugate pair for fc € (8 — 6-\/2,8 -I- 
6-\/2). For our calculations, we have taken k = 2 which 
shows that at e* = = 2, Ai crosses the zero line 

as can be seen in Fig. [3] a) (red symbols). For higher 
values of £ > stationary states {x*+,y*_^_,u*+) (cir¬ 
cles) in Fig. [3] (a) for e > £*(= Etc) get stabilized via 
a transcritical bifurcation where {x*,y*) = (1,0) and 
(x*_|_,y*_|_,u*+) exchange their stability. This again is 
quite unexpected since the feedback is designed to stabi¬ 
lize {x*, y*) = (1, 0) for higher £ values. A brief discussion 
regarding the behavior of this system in the {e, k) plane 
is available in Appendix ICl 

For the origin at {x*,y*) = (0,0), the bifurcation di¬ 
agram (black dots) is shown in Fig. |3]b). Appropriate 


transient trajectories corresponding to different augmen¬ 
tation regimes are also shown in related Figs. |3] (b.l), 
(b.2). We observe bistability for a range of lower £ values 
before the origin gets stabilized. The stationary solutions 
obtained in the bistable regime are given by 


x°±=±\ 1- 


y ± = X ±-X ± 

.yl± 

e 


k — 

3 


u ± = 


(16) 


Transient trajectories in this regime demonstrating 
the two observed stationary solutions are shown in 
Fig. [3] (b.l). These solutions approach and collapse at 

the origin at a pitchfork bifurcation for epp = 

(=1 for A: = 2 in this case) beyond which the solutions 
{x°±Ty°± tU°±) become imaginary and the origin is the 
only stable real stationary solution. A transient trajec¬ 
tory in this parameter regime is shown in Fig. [3] (b.2). 
The characteristic equation for the origin is, 

A^ + kX^ + X{2e^ -1) + 2£^ - k = Q, (17) 



which has all negative eigenvalues for 2£^ — /c > 0 giv¬ 
ing us a stability regime of £ > \fkj2 and a transition 
value oi Epp = e* = \fkj2 = 1 (since k = 2) when the 
eigenvalue/s cross the zero axis. From Eq. 1171 we get 
Ai = —1 in the large £ limit. It is numerically observed 
here that the discriminant A < 0 in this range and there¬ 
fore Re(A 2 , 3 ) = (1 — A:)/2 = —0.5 and consequently, the 
origin is stable in the large £ limit. The largest eigenvalue 
for the origin is plotted as red symbols in Fig.[3](b) which 
shows the changes in the stability of the origin from un¬ 
stable in £ G (0,1) to stable V £ > 1. For a discussion 
regarding the system behavior in the (£, k) plane, please 
see Appendix ICl 

For {x*,y*) = (—1,0), the bifurcation diagram (black 
dots) is shown in Fig.[3]c). Appropriate transient trajec¬ 
tories corresponding to different augmentation regimes 
are also shown in related Figs. [3] (c.l), (c.2), (c.3). Since 
this solution is a symmetric counterpart of {x*,y*) = 
(1,0), the corresponding analysis similarly carries over 
in this case. 

These simple examples demonstrate the fact that tar¬ 
geting the required stationary solutions using linear aug¬ 
mentation is not quite straightforward and the procedure 
is quite sensitive to how the systems are augmented, the 
stationary solutions being targeted and to the properties 
of systems. In the following, results for a specific class 
of dissipative dynamical systems are presented to further 
highlight these limitations. 


C. Dissipative predator—prey models 

Considering predator-prey population models, general 
evolution equations for these systems with logistic prey 
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FIG. 4. (Color online) Bifurcation diagram (black dots), largest eigenvalues (red circles) and time series (insets for specific 
e values) for predator-prey systems with H II (left column) and H III (right column) functional responses. Top row: For 
augmented prey, insets show the time series of x,y and u for two different e values before (with oscillatory u) and after 
the stabilization of (x* ,y*) = (0,0) (with u* =0). The Hopf bifurcation is marked as H, and the transcritical bifurcation 
points have been highlighted by Ti, and T 2 respectively. Middle row: For augmented predator, the systems exhibit oscillatory 
behavior similar to augmented preys (not shown) for low e before they settle on the stationary solution {x*,y*) = [K,Q) 
with u* = 0; where the preys reach their carrying capacity in the absence of predators for higher e values. Bottom row: For 
augmented predator and prey, for low e, the systems are oscillatory (not shown). With increasing e, both the systems lose the 
oscillatory behavior and all trajectories escape to infinity (gap in the bifurcation diagram with no black dots). For higher e 
values, unrealistic stationary solutions where either the preys exceed their carrying capacity (a;* > K) with negative predator 
populations (y* < 0) (H II and H III), or where the prey populations are negative with small positive predator population (for 
H II) and u* yf 0 get stabilized. 


growth can be written as, 

X = rx{l - x/K) - f{x)y, 

V = {pf{x) - 7 )y, (18) 

where x and y correspond to prey and predator popula¬ 
tions respectively and the parameters r, K, p, and 7 are 
positive. Considering the evolution equation for preys, 
the first term rx{l — x/K) represents the logistic growth 
rate of the prey species with the maximum growth rate of 
r and carrying capacity K which is the maximum popu¬ 
lation size that the environment can sustain indefinitely. 
The second term f{x)y corresponds to the prey mortal¬ 
ity via predation. f{x) is the functional response gov¬ 
erning the rate of per capita prey consumption by the 
predators The parameter p governs the biomass 

conversion efficiency for the predators in the sense of 
how many predators are added to the population via 


predation, and 7 is the intrinsic predator mortality pa¬ 
rameter. One of the stationary solutions of this system 
corresponds to vanishing predator-prey populations, i.e. 
{x*,y*) = (0,0). The other stationary solutions are de¬ 
pendent on the type of functional response considered. 
Most commonly employed f{x) forms in such models are 
the Holling type with the following general expressions: 


1. f{x) = ax for Holling type I response which is 
identical to the predation in the Lotka-Volterra 


case 


2. fix) = 


for Holling type H (Michaelis- 


ib + x) 

Menten kinetics), using which, Eq. [TSl gives the 
Rosenzweig-MacArthur model , 


3. fix) = 


ax 


(&2 


for Holling type HI (Hill equa- 
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tion type), usin g wh ich, Eg. [TSl gives the Truscott- 
Brindley model |40l | which is used in modeling phy¬ 
toplankton and zooplankton interactions leading to 
harmful algal blooms, 

and consequently, corresponding stationary solutions can 
be obtained. The parameter a in expressions above cor¬ 
responds to the maximum per capita predation rate, 
and b is the half saturation constant governing how 
quickly the predators attain their maximum consump¬ 
tion rate. In the following, we will have a closer look 
at the stability properties of the trivial stationary so¬ 
lution {x*,y*) = (0,0): origin. Considering a general 
augmented population model, 

X = rx{l — x/K) — f{x)y + £iu, 

y = (pf{x) - 7 )y -h S 2 M, 

ii =-ku-ei{x - X*) - e 2 {y - y*), (19) 

it turns out that the Jacobian for this system is identi¬ 
cal for all three functional responses at the origin. The 
identical characteristic equation therefore is, 

(r - A )(7 -I- X){k -I- A) -I- e 2 ^(r - A) - £ 1^(7 -I- A) = 0. 

( 20 ) 

For £1 = £2 = 0, we obtain the eigenvalues as Ai = 
r, A 2 = — 7 , and A 3 = —k where Ai and A 2 are the eigen¬ 
values for the original system in Eq. [TH] implying that 
the origin is unstable, and A 3 corresponds to the expo¬ 
nentially decaying augmentation variable u. Since the 
Holling type I case with insatiable predators is quite un¬ 
realistic, we will focus here on systems with Holling type 
II (H II) and III (H III) behaviors. In the following anal¬ 
ysis, the parameter values are fixed at: r = 0.5, K = 0.5, 
a = 1/3, b = 1/15, p = 0.5, 7 = 0.1 for the H II [U 
system, and r = 0.43, K = 1, a = 1, 6 = 0.053, p = 0.05, 
7 = 0.028 for H III (^,[1^. Let us now look at different 
augmentation situations: 

For £1 = £ and £2 = 0, i.e. only prey populations are 
augmented, substituting these values in Eg. [20l gives. 

(7 + A)[(r-A)(A: + A)-£2]=0. (21) 

Since one of the roots A = —7 is independent of £ there¬ 
fore the remaining roots of this equation determine the 
stability of the origin. The remaining two roots are 
A± = — (/c — r)/2 ± \J(k + r)2 — 4£2/2 out of which, 
A_ < 0, V £. It is easily verifiable that the eigenvalue 
A+ (which also is the largest) is positive V £ < and 
crosses the zero axis at e* = leading to all nega¬ 
tive eigenvalues and hence a stable origin. This is quite 
similar to the harmonic oscillator case where increas¬ 
ing/decreasing the value of the decay parameter k could 
increase/decrease the threshold value of stable —?► unsta¬ 
ble transition (unstable —?► stable in this case). Further¬ 
more, in the large e limit, we obtain the largest eigenvalue 
Ai = —7 and therefore the origin is stable in this reg ime. 
Fig. m top row shows the bifurcation diagram and 
the largest eigenvalue behavior for H II (left) and III 


(right). The unstable —>■ stable transition in both these 
systems with increasing coupling values can be seen in 
the figure. Although for the H III system in the regime 
£ > 0.927, for r = 0.43, k = 2), certain initial con¬ 

ditions lead to the trajectories escaping to infinity (not 
shown) which accounts for the missing dots in the bi¬ 
furcation figure. Since x and y are population variables 
by definition, population models are constrained to work 
for non-negative values of x and y respectively. What 
we observe here is that the augmentation forces the prey 
populations into taking negative values which leads to a 
breakdown in the model constraints and the logistic func¬ 
tion in the rate equation of x leads to diverging solutions 
as time increases. For H II system this appears not to 
be the case and all considered initial conditions lead to a 
stable origin V £ > \/^(= I, for r = 0.5, k = 2). 

For £ 1=0 and £2 = £, i.e. only the predator popula¬ 
tions are augmented, substituting these values in Eq. [201 
give, 

(r — A )[(7 -|- X){k -|- A) -I- £^] = 0, (22) 


and we see that an eigenvalue A = r is always positive 
since r > 0 , and therefore this setup will never stabilize 
the origin. The remaining eigenvalues are A± = — (7 + 
k)/2 ± \/(k — j) + 4e^f2. In Fig. |4j second row, for low 
£ values, the systems exhibit periodic behavior similar to 
the one shown for the augmented prey case. For higher 
values of £, both H II and H III settle on a stationary 
solution of the original system (x* ,y*) = (K, 0) which we 
did not intent to stabilize. For this solution, the preys 
exist at their carrying capacity and the predators vanish. 
For H II and H III, the carrying capacities considered for 
simulations are K = 0.5 and K = 1 respectively, and 
hence the observations in Fig. |4] (middle row). 

For £1 = £2 = £, i.e. both prey and predator popula¬ 
tions are augmented, substituting these values in Eq. 
and rearranging terms gives. 


X^ + (^ + k — r)A^ -I- (2£^ + (7 — r)k — rj)X 

+£^(j — r) — rjk = 0. (23) 


Using the RHC [S^, one of the conditions for this equa- 

I 7^^ /c 

tion to have all negative roots is £ > W 7 - 7 which 

\ il-r) ^ 

is impossible to achieve since r > 7 . Therefore, this 
setup will not stabilize the origin either. Fig. 2] bottom 
row shows the behavior of H II and H III. For smaller 
£ values, these systems exhibit periodic behavior similar 
to the augmented prey. On increasing e further, systems 
enter a regime where all considered initial conditions lead 
to escaping trajectories. The reason behind this behavior 
here again is due to a breakdown in modeling constraints. 
Examination of transient trajectories reveals that aug¬ 
mentation in this case is forcing the predator populations 
into y < 0 axis which leads to a breakdown in the model, 
thereby initiating a positive feedback loop in the prey 
populations leading to the diverging behaviors observed 
in simulations. Beyond this regime for higher values of 






HII H III 



e e 

FIG. 5. (Color online) Different dynamical regimes for H II 
and H III systems are marked as A, B and C in the e, k plane. 
A is the regime of periodic dynamics, in B stationary solutions 
of the coupled system are stable and C is the regime of stable 
origin. The boundaries between A—>-B and B—>-0 are the loci 
of the reverse Hopf bifurcation H and the second transcritical 
bifurcation T 2 respectively (as in Fig. [4] top row). 


e, H II system exhibits bistability between different sta¬ 
tionary solutions where in one case, preys exceed their 
carrying capacity {x* > K) and the predator popula¬ 
tions are negative {y* < 0), and in the other case, the 
prey populations are negative (x* < 0) and predators as¬ 
sume a small positive value. It is important to note yet 
again that these solutions are impractical because the 
populations cannot exist above their carrying capacities 
nor can they take negative values under realistic model¬ 
ing constraints. For H III system, we only observe the 
equilibrium solutions with x* > K and y* < 0 (see inset). 
In both the cases, we have a non vanishing u* > 0 and 
therefore these solutions exist due to augmentation and 
cannot be observed otherwise. Following this analysis, 
we can conclude that augmenting the prey is the correct 
strategy to stabilize of the origin and the other coupling 
schemes can lead to complicated dynamics. Even though 
the analysis here is limited to the origin, we can expect 
these behaviors to be quite general with respect to other 
stationary solutions as well. 

Now considering applications, as already mentioned, 
origin corresponds to an equilibrium for which the preda¬ 
tors and the preys vanish. Persistence of populations 
for a proper ecosystem function is very imperative and 
has been studied extensively from several perspectives, 
contributing towards a better understanding of the pro¬ 
cesses leading to species extinction Knowledge 

regarding these processes can help in devising procedures 
which can contribute towards better species conservation 
efforts. For the simple models considered in the previous 
analysis, it is clear that either by coupling the system 
appropriately or by using specific parameter values for 
k and e, we can avoid stabilizing the origin. For in¬ 
stance, considering the prey augmented case, for low e, 
the systems exhibit periodic oscillations. On increasing 
the coupling strength, stationary solutions of the aug¬ 


mented system {x* > 0,y* > 0,u* < 0), satisfying 

rx*(l - x*/K) - f{x*)y* + em* = 0, 
ipfix*) - l)y* = 0 , 

-ku* - ex* = 0, (24) 

get stabilized through a reverse Hopf bifurcation (marked 
as H in Figs.|4](top row)). For these stationary solutions, 
the value of x* stays constant while y* and u* = —£x*/k 
show a variation for a range of e values (plateau between 
H and Ti in Fig. 0] (top row)). It is also important to 
note that some initial conditions in this regime can lead 
to trajectories escaping to infinity. This branch of so¬ 
lutions undergoes a transcritical bifurcation (marked Ti 
in Fig. m (top row)) where it exchanges stability with 
another branch of solutions with u* —>■ 0 for increasing 
e. At £* = \frk, M* = 0 and the predator-prey system 
effectively decouples from augmentation which is accom¬ 
panied by another transcritical bifurcation {T 2 in Fig. 0] 
(top row)) between the continuing branch of stationary 
solutions (x* > 0, y* > 0 ,m* —>■ 0) and the origin. In 
£ > £* regime, origin is the only dynamical attractor. 
Fig. [ 5 ] shows the parameter scans for H II and H III sys¬ 
tems highlighting these different dynamical regimes. In 
region A these systems exhibit periodic behavior and the 
boundary between A and B is the locus of the Hopf bifur¬ 
cation in the e, k plane, which leads to the stabilization 
of stationary solutions {x* > 0,y* > 0,u* < 0). B cor¬ 
responds to the regime where stable stationary solutions 
{x* > 0,y* > 0,u* < 0) and {x* > 0 , 2 /* > 0,u* 0) are 

observed and the boundary between B and C is the lo¬ 
cus of the second transcritical bifurcation T 2 which leads 
to the stabilization of the origin. Therefore by using ap¬ 
propriate values of e and k, we can keep the system in 
either a periodic state, or a stationary state with non 
vanishing populations and can expect this procedure to 
work in experiments and be robust with respect to noise; 
Ref. experimentally stabilized a stationary solution 
in an electronic Lorenz system at permitted noise level. 
Furthermore, in the other instances of augmented preda¬ 
tors, or augmented predators and preys we already ob¬ 
serve a complete lack of origin stabilization. Therefore, 
we can employ these schemes as well to avoid stabiliz¬ 
ing the origin but one needs to be careful since these 
cases can lead to other complications as discussed. An¬ 
other useful application for these observations could be 
in cases where maximization of prey yield is required. 
Augmenting the predator populations is seen to stabi¬ 
lize the equilibrium where the prey populations exist at 
their carrying capacity and the pre dators vanish. This 
can find ^plications in fisheries |48l . , algae fuel gen¬ 

eration [^ : where maximal sustainable yields are 
crucial, and also in biomedical research, for e.g. in HIV- 
1 infection models [52| where a portion of human immune 
system i.e. activated CD4+ T cells are the primary tar¬ 
get of the HIV-1 infection which can be modeled 

via predator-prey dynamics. 
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IV. SUMMARY 

In this work we studied the general ability of linear 
augmentation towards stabilizing desired stationary 
solutions of oscillatory systems. Through some simple 
examples discussed in this paper, it is clear that the 
effectiveness of this scheme is quite sensitive to the aug¬ 
mentation parameters, the class of oscillatory systems 
considered, the stationary solutions to stabilize and 
also on the way the systems are augmented. Therefore, 
although the simplicity of linear augmentation makes 
it a very compelling choice for applications, a careful 
analysis is required to test the system for potential 
pitfalls associated with the scheme. As highlighted by 
the examples, apart from failing to target the appropri¬ 
ate stationary solutions, linear augmentation can also 
lead to other complicated dynamical situations which 
include escaping trajectories, stabilization of unintended 
stationary solutions or the stabilization of stationary 
solutions which are not permitted under the modeling 
constraints; preys existing above their carrying capaci¬ 
ties and negative predator populations in Sec. IIII Cl for 
instance. Nevertheless, one can find ways to exploit the 
failures of the scheme in applications. Although we can 
expect to see these results in experiments, an in-depth 
study of this procedure in presence of noise, and also for 
larger systems is required. Extending on the results in 
the ecological context, one needs to check the process 
behavior in presence of multiple preys and predators, for 
a food chain, and also for other functional responses [55l| . 
Furthermore, linear augmentation has been proposed as 
a mechanism to control bistability but how it fares 
in managing more general instances of multistability 
including extreme multistability [H, is still an open 
question and will be addressed in subsequent studies (55l| . 
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Appendix A: Converging and diverging trajectories 
in augmented harmonic oscillator 



FIG. 6. (Color online) Bifurcation diagram (black dots) along 
with the largest eigenvalue (red symbols) as a function of £ in 


the top row. e* = lo 


OJ^-1 


(~ 1.63) marks the coupling be¬ 


yond which all initial conditions lead to escaping trajectories. 
Transient trajectories shown for e{= 0.5) < £* (left bottom) 
and e(= 1.75) > e* (right bottom). 


\/efei-4(e? + u}'^)/2. For partially augmented cases 
a = 0 and /3 = -\/^^4(e^"+w^/2 or P = ioj for £2 = 
0 , £1 ^ 0 and £1 = 0 , £2 7 ^ 0 respectively. 

For identical augmentation £1 = £2 = £, we get ID = 
-|- e'^D + + w^) and Eg. I All reads 

Da; = U{e, k, t)[= £{1 -I- £^ — k}u{t)]. (A2) 

The roots of the auxilliary equation in this case are m± = 
a±P with a = —£^/2 and P = \J— 4(£^ -|- a;^)/2. For 
imaginary /3, the transient solution for Eq. IA2I can be 
expressed as. 


Xg (t) = Axi {t) + Bx 2 {t ), (A3) 


which is independent of the forcing term U{e,k,t) with 
xi(t) = exp (at) cos pt, and X 2 {t) = exp (at) sinpt. Con¬ 
sequently, the steady state solution can be obtained by 
using the Laplace and inverse Laplace transformations 
giving. 


1 

Xst{t) =sin.{VL{t — x))U{e,k,x)dx, 
“ Jo 


(A4) 


Augmented harmonic oscillator dynamics from Eq. [S] 
can also be expressed in form of a second order ODE as. 

Da; = [/(£i,£ 2 ,fc,t)[= {£2 -I-£ 2 ^ 1 ^ -£ifc}u(t)], (Al) 

where the derivative operator D = D^-|-£i£ 2 D-l-(£^-|-w^) 

dJ 

with D® = —r, i = 1,2 in this case. This equa- 

dP 

tion corresponds to a driven harmonic oscillator with 
frequency (ef -I- w^) and a damping coefficient £i£ 2 . 
The roots of the auxiliary equation for the operator 
D are m± = a ± P where a = —£i£ 2/2 and P = 


where fl = \/oJq — 7 ^, = £^ -b and 7 = £^/ 2 . 

Now at this point, we do not know the exact expres¬ 
sion for U{e,k,t)- Considering the transient behavior 
of trajectories in partially/fully augmented system, we 
clearly observe that they possess an exponentially de¬ 
caying/diverging envelop (see Fig. [T] (inset) and Fig. [5] 
bottom row). Based on these observations, assuming 
U{e, k, t) = ao exp (kmt) where both ao,km are functions 
of £ and k, and solving Eg. I A4I gives 

= (A5) 
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FIG. 7. (Color online) Largest eigenvalues and estimated val¬ 
ues of km as functions of e, (a) for the system augmented in 
X and (b) for the fully augmented system. 



FIG. 8. (Color online) Maximal Lyapunov exponent for in¬ 
creasing and decreasing values of e for the fully augmented 
DufHng system with (x* = l,y* =0) demonstrating hystere¬ 
sis. Calculations for increasing and decreasing e are marked 
by red and black arrows respectively. 

From this expression we see that the trajectories will ex¬ 
ponentially decay to the origin V km < 0 and diverge V 
km ^ 0 - 

Fig.[7]shows the numerical estimation of km along with 
the largest eigenvalue of the Jacobian Ai at the origin as 
a function of e; for partially (Fig. [7] (a)) and fully aug¬ 
mented cases (Fig. [7] (b)). km here was calculated as the 
average rate of convergence/divergence in the Euclidean 
distance of the current systems’ state from its previous 
state, for every time step along the trajectory. These re¬ 
sults suggest that km = and this observation has some 
interesting consequences. For km = < 0, we have a 

case of a harmonic oscillator being driven by an exponen¬ 
tially decaying force and consequently the system settles 
on the origin as t > oo. Similarly the other case of an 
unstable origin [km = Ai > 0) is equivalent to the os¬ 
cillator under the influence of an exponentially diverging 
force (energy being pumped into the system) which leads 
to diverging trajectories as time increases. 


Appendix B: Hysteresis in augmented Duffing model 

For the bistable fully identically augmented Duffing 
system in Eq. [8] with x* = l,y* = 0, the largest Lya¬ 
punov exponent was calculated for increasing and de¬ 
creasing values of augmentation strength e. Starting with 
initial conditions leading to the solutions {x*-,y*z*-) 



0 2 4 6 0 1 2 3 

e e 


FfC. 9. (Color online) Behavior of fully augmented Duffing 
system in the (e, k) plane. Parameter regimes marked in grey 
(A) correspond to a successful stabilization of the intended 
stationary solution; (1,0) in (a) and (0,0) in (b). Regimes 
of bistability are marked with blue dots and B corresponds 
to parameter values for which other stationary solutions are 
stable. 

at e = 0.1, the initial conditions for the next calculation 
at e = O.l-bfe were taken as the final values of x, y, u from 
the previous calculation for e = 0.1, with 6e = 0.001 and 
so on for the entire range in the forward direction. Simi¬ 
larly for backwards calculation, the process was repeated 
starting from e = 0.7 where (1,0) is the only stable so¬ 
lution with Se = —0.001. The results of the calculation 
are shown in Fig. [5] and as one would expect, this system 
exhibits hysteresis in the interval of bistablity. 


Appendix C: Fully augmented Duffing system: 
behavior in {e, k) plane 

Different dynamical regimes for the Duffing system in 
Eq.[8]with {x*,y*) = (1,0) and (0,0) are shown in Fig.|9l 
The grey areas marked as A correspond to the regimes 
where the intended stationary solutions are successfully 
stabilized, namely (1,0) and (0,0) in Figs. IHKa) and (b) 
respectively. 

For X* = l,y* = 0 in Fig. Ela): blue dots for lower e 
values highlight the regimes of bistability where solutions 
{x*-,y*_,z*-) (from Eq. (TS]) and (1,0) coexist. Solu¬ 
tions (a:*_, t/*_, z*_) vanish via a saddle node bifurca¬ 
tion after colliding with the unstable branch of solutions 
{x*j^,y*j^,z*+) (again from Eq.[T31) and the boundary of 
the blue dot regime gives the locus of this saddle node 
bifurcation; parabolic function Fsn{£, k){= 5e^ — fc) = 0, 
estimated from the expression for x* ± from Eq. [T3] as 
the limiting value of k and e to get real solutions, ob¬ 
tained by equating the discriminant in the expression 
of X* ± to zero. The boundary between regimes A and 
B corresponds to the locus of the transcritical bifur¬ 
cation between (1,0) and [x* + ,y*z*+) given by the 
zero crossing of the largest eigenvalue for (1,0) which 
satisfies Ftc{£, k){= 2e^ — k) = 0. In region B either 
the now stable coupling dependent stationary solutions 
{x* + ,y*z*+), or escaping trajectories are observed. 

Similarly for x* = 0,y* = 0 in Fig. HKb): B high¬ 
lights the regime where the system settles on the so¬ 
lutions (x°j^ , u°+). The blue dots correspond to 

the initial conditions leading to stationary solutions 
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{x°-,y°_,u°-). Since the system is bistable in this 
regime, we should expect the entire region B to be 
densely filled with these blue dots but that is not the 
case. The reason behind this behavior is a difference in 
the relative basin size of these two solutions; number of 
initial conditions leading to (a:°+, y°^,u°+) is more than 
the ones which lead to {x°-,y°_,u°-). This difference is 


even more pronounced for higher values of k. Both these 
solutions vanish via a pitchfork bifurcation and the lo¬ 
cus of this bifurcation which separates regimes B and A 
can be traced by the function Fpf{s, k){= 2e^ — k) = 0 
which is obtained by equating the discriminant in the 
expression of x°± from Eq. [TClto zero. 
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